Numerical Investigation of the Dynamics of a Thin 
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The equilibrium dynamics of a thin film type II superconductor with spherical geometry are 
investigated numerically in a simulation based on the lowest Landau level approximation to the 
time-dependent Ginzburg-Landau equation. Both the static and time-dependent density-density 
correlation functions of the superconducting order parameter have been investigated for systems 
with varying amounts of quenched random disorder. As the temperature is lowered it is found that 
the correlation length, the length-scale over which the vortices have short-range crystalline order, 
increases but the introduction of quenched random disorder reduces this correlation length. We see 
no signs of a phase transition in either the pure or the disordered case. For the disordered system 
there is no evidence for the existence of a Bragg glass phase with quasi long-range correlations. The 
dynamics in both the pure and disordered systems is activated, and the barrier of the relaxation 
mechanism grows linearly with the correlation length. The self-diffusion time scale of the vortices was 
also measured and has the same temperature dependence as that of the longest time scales found 
in the time dependent density-density correlation function. The dominant relaxation mechanism 
observed is a change in orientation of a correlated region of size of the correlation length. A scaling 
argument is given to explain the value of the barrier exponent. 
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I. INTRODUCTION 

The nature of the mixed state of thin films of type 
II superconductors is despite considerable experimental, 
theoretical and numerical research still unclear. In the 
mean- field limit, which describes the ground state of the 
system, the vortices are known to form into a triangular 
Abrikosov lattice jlj]. Also it is generally accepted that 
close to the mean-field transition line to the normal state, 
the H C 2 line, thermal fluctuations destroy that hexago- 
nal order and lead to a vortex liquid state B. However, 
the phase diagram between the mean-field H c x and H C 2 
lines remains unclear. The Kostcrlitz-Thoulcss-Halperin- 
Nelson (KTHN) theory of a two dimensional (2D) melt- 
ing U has been applied to superconductors M M yield- 
ing a phase diagram with a solid, a hexatic and a liq- 
uid phase. Other theoretical work raises a doubt as to 
whether there is any phase transition to the mixed state 
and suggests that a flux lattice exists only in the zero 
temperature limit and that the liquid phase exists at any 
non-zero temperature |J Q. It has long been predicted 
that in the presence of disorder any long range crystalline 
order of a vortex lattice phase is destroyed in less than 
four dimensions ||. In the glass like state which is ex- 
pected to form instead, the length scale of the short-range 
crystalline order has been predicted by Ovchinnikov and 
Larkin jj| for the zero-temperature (mean-field) limit. 
More recently the existence of a Bragg glass phase with 
quasi long range order has been suggested ]Tc| . The non- 
perturbative analytical method used by Yeo and Moore 
0, which yields a 2D flux liquid at all temperatures, has 
also been applied to the case of quenched random dis- 



order, which is found to just reduce the extent of short- 
range crystalline order present in the liquid state. We 
shall find good agreement of our numerical results with 
this work, but none with the KTHN picture, or in the 
presence of disorder with the Bragg glass scenario. 

A lot of numerical work has been done on clean sys- 
tems, some indicating a first order 2D vortex melting 
transition |ll[] and some failing to see any kind of 
phase transition |l3fl fljj ] ■ The experimental evidence 
on superconducting films is also contradictory. Whilst a 
first order 2D melting transition has never been observed, 
current-voltage measurements have provided some evi- 
dence for second order 2D melting |l6| as well as failing 
to provide evidence for any transition in other experi- 
ments |lj]] . Yazdani et al. |l8j have reported evidence for 
2D melting from measurements of ac penetration depth 
in thin films of a- MoGe. Theunissen et al. find indi- 
cation of a melting as well as a hexatic to liquid transition 
in measurements of vortex viscosity in NbGe. 

Most of this experimental evidence is on transport phe- 
nomena in samples that have at least weak disorder or 
surface pinning due to crystal defects or impurities. This 
focuses interest on the dynamics of the mixed state and 
the influence of quenched random disorder. For the de- 
scription of ac transport phenomena equilibrium dynam- 
ics is relevant. It is also vital for the whole picture of 
the mixed state. It has been suggested that the apparent 
second order freezing transition deduced from changes 
in vortex liquid viscosity or sample resistance may in- 
stead be due to an exponentially fast increase of the 
relaxation time scales at low temperatures. p0[ . This 
behavior, which our work confirms, indicates activated 
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dynamics, in which the relaxational time scales increase 
exponentially with the correlation length. 

In this paper we present a numerical investigation 
into relaxational time scales of density correlations in 
the vortex liquid and diffusion time scales of the vor- 
tices in a thin film. A Langevin dynamics simulation is 
used with a phenomenological Ginzburg-Landau Hamil- 
tonian. To simulate quenched random disorder, we 
add a Gaussianly-distributed random contribution to the 
mean-field transition temperature in the Hamiltonian. 
The equilibrium dynamics of the superconductor is in- 
vestigated by measurements of the time scales of density 
correlation decay and vortex self-diffusion. We find that 
the time scales of the two different processes both show 
the same exponential growth with the translational cor- 
relation length in the system and are due to the same 
activated relaxation process. 

The paper is organized as follows. In Section || wc 
shall give a short summary of the background of the 
phenomenological Ginzburg-Landau theory on which our 
work is based. We then give a description of how we solve 
the time dependent Ginzburg-Landau equation numeri- 
cally in Section III. Section IV reports our work on the 
density correlations and their relaxation times, followed 
by Section ^ which reports on self-diffusion of the vor- 
tices. In Section VI we put together our results from 
the previous sections and interpret them. We also iden- 



tify the underlying dynamic processes. In Section VII wc 
summarize and discuss our results. 



II. THEORETICAL FRAMEWORK 

Our simulation is based on the phenomenological 
Ginzburg-Landau theory in the approximation of a uni- 
form magnetic field B. The free energy functional for a 
clean sample is given by 

F = J d\ (aH 2 + + 2^|(-»ftV - 2eA)Vf) , 

(1) 

where ip is the order parameter representing the macro- 
scopic wave function of the superconducting electrons 
and B = V x A. In first approximation a(T) = a'(T — 
T c ) and (3(T) is a constant, a', > 0. 

For the case of quenched disorder a random local po- 
tential is added to the free energy: 

F dis = J d 3 r9(r)|V>(r)| 2 , (2) 

where 0(r) is real and Gaussian distributed with 

(6(r)) = 0, (3) 



Angular brackets denote thermal averages and 8 is the 
three dimensional Dirac delta function. A is the mea- 
sure of the strength of the disorder. 

The simulation follows Langevin dynamics, described 
by the time dependent Ginzburg-Landau equation. 



ay»(r,t) 

dt 



dF 



(•5) 



(e(r)e(r')) = AJ(r-r') 



(4) 



where F is defined in Eq. (|l|) and £ is Gaussian white 
noise of strength 2rfcsT, i.e. 

(e(r,t)^r , ,t , )) = 2Tk B TS(r-r')5(t-t') . 



III. NUMERICAL SOLUTION OF THE TIME 
DEPENDENT GINZBURG-LANDAU EQUATION 

We base our simulation on a thin film superconductor 
with the following experimentally realizable two dimen- 
sional geometry. A thin film superconducting sphere of 
radius R and thickness d is placed in a radial magnetic 
field of magnitude B which emerges from an infinitely 
long, thin solenoid whose end is at the center of the 
sphere. This system has been investigated numerically 
before, using Monte Carlo dynamics Q Q Q. The 
reasons for our preference of this geometry to the more 
widely used geometry of a plane with periodic boundary 
conditions have been presented in detail by Dodgson and 
Moore jllj] . The main advantage is that the spherical ge- 
ometry guarantees full translational symmetry, which pe- 
riodic boundary conditions do not, while having the same 
thermodynamic limit. Periodic boundary conditions im- 
pose an artificial pinning potential on the vortices which 
makes them unsuitable for investigating transport and 
dynamical phenomena. 

We assume that in a strong type II superconductor 
where the magnetic penetration depth is larger than the 
coherence length of the superconducting wave function 
by orders of magnitude, fluctuations in the magnetic in- 
duction are negligible. Flux quantization requires that 
the B field at the sphere be such that An R 2 B = N$ , 
where $o = h/2e is the flux quantum and N is the num- 
ber of vortices. This fixes R = y/N/2 with the unit of 
length being l m = \J <Pq/2ttB, which is proportional to 
the nearest neighbor distance of vortices in the Abrikosov 
lattice. 

We shall make the usual approximation of retaining 
only the lowest Landau level (LLL). This is the approx- 
imation that Abrikosov used to first describe the vortex 
lattice state at mean-field level M. Defining D as the 
gauge invariant momentum operator, D = —iKV — 2eA, 
the LLL is the eigenspace of D 2 associated with its lowest 
eigenvalue 2eBh. The LLL approximation is tradition- 
ally believed to hold near the upper critical field H c2 . 
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However, due to renormalization effects, it actually de- 



scribes a large portion of the vortex liquid regime |21|. 

In spherical geometry an orthonormal basis of the LLL 
is || 

iMM) = (4vri? 2 )- 1 / 2 A m e lm *sin m (0/2)cos Ar -™(0/2), 

(6) 

where N is the number of vortices, < to < N, and 
A m = B(m + l,N - m + l) -1 / 2 . B stands for the Beta 
function, given for natural numbers n,m by 



Bin, to) 



(n - 1)! (to - 1)! 
[n + rn — 1)! 



The order parameter in the LLL approximation can be 
expanded in the above basis set: 



N 



ip(9, 4>) = Q } v m tp m (0, 4>), 



m— 



(7) 



where Q = (<$> k B T / fidB) 1 ^ Q. 

The free energy Hamiltonian for our system is given 
as a sum of the Hamiltonian of the clean system and the 
contribution from disorder, Ti = Ti c i + Tidis, where Ti c i 
and Tidis are the free energy terms from Eqs. (Q) and (|^) 
respectively. If the order parameter is now expanded in 
the LLL eigenstates, we can carry out the spatial integral 
to express the Hamiltonian of a system without disorder 
in terms of the LLL coefficients v m : 



v rr w N 1 2N 

m—0 p=0 

The effective temperature parameter is given by 



(XT 



dQ 2 



a(T) 



eBh 



(8) 



(9) 



Note that ar — corresponds to the mean- field H C 2(T) 
line and ax = — oo to T = 0. In the quartic term, 



N 
9=0 



6(p ~ q)8{N - (p - q)) fN(q,P~ q) v q v p - q , 



(10) 

where 9 is the Heaviside step function and 

f N {m,n) = (A m A n B(m + n + l, 2N - to - n + 1)) 1/2 . 

To express the disorder contribution Tidis in terms of 
the LLL coefficients v m , we expand the random Gaussian 
disorder on the sphere in normalized spherical harmonics 

ym. 



7, rp OO I 



with a[™* 



To satisfy Eqs. 



© and (|), the 
real and imaginary parts of the a™ are drawn inde- 
pendently from a Gaussian distribution with variance 
D = (dQ 2 /k B T) 2 A/2 for m^O and the a from a Gaus- 
sian with variance 2D. Then the part of the Hamiltonian 
due to disorder can then be expressed 



Ti d is{{v m }) _ V o, v*v T lp ~ ql 



knT 



0,9=0 l=\p-q\ 



where 



ir n = / d'rYrr n+m ^n 



(p+q-\p-q\)/2 

(11) 

(12) 



y 4tt(1 — to)! to! 
m—l, to+Z+1, n+TO+1, 



x 3-F2 



to+1, N+m+2 ; 1 



with 3F2 a generalized hypergeometric function. This 
term of the Hamiltonian limits our simulation in the case 
of disorder to a relatively small system size, as it requires 
finding and storing ~ N 3 /6 values of 3^2- However, as 
the correlation length is very much reduced by disorder, 
the effects of the finite system size are for the same Oct 
much reduced compared with the non-disordered case. 

The time dependent Ginzburg-Landau equation, dis- 
cretized in time, which drives our simulation is now: 



v m (t + At) - v m (t) 

^d{H c i{t)+H dls (t)) 



(13) 



-Atr- 



dvl 



AtU{t) 



1=0 m=-l 



with H c i(t) and Hdis(t) as defined in Eqs. (g) and (p. l|). 

This equation ensures that thermal averages can be 
replaced by time averages over successive measurements 
from the simulation. It leads to the correct canonical 
distribution after infinite time in the limit At — > 
if real and imaginary parts of the £ m are drawn inde- 
pendently from a Gaussian probability distribution with 
variance a 2 /At with a 2 = 2Tk B T @. However, the 
finite time steps always lead to a certain small deviation 
from the correct distribution, which does not arise in a 
Monte Carlo algorithm H]. We have tested our sim- 
ulation's static thermal averages for varying quantities 
like energy and entropy against Monte Carlo averages 
from the code used in Ref. Jl5| , and chose our time steps 
At = 0.15, which reduces deviations to less than 1.5%. 

The quenched disorder is self-averaging for the case of 
an infinite system. In a finite system though, the par- 
ticular random choice of disorder will obviously influence 
the results. To reduce this effect we run the simulation 
between 10 and 20 times with different sets of random 
aY 1 , and average over these different measurements. 
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IV. THE DENSITY-DENSITY CORRELATOR 



C(k,0) 



To investigate the extent of order and the time scales 
of fluctuations in the system, we compute the connected 
part of the density-density correlator. This correlator 
measures the correlations of the magnitude of the order 
parameter in space and time. The density-density cor- 
relator carries essentially the same information as the 
translational correlation function of the vortex positions 
jL3[. However, it is far easier to compute than the latter, 
because it does not involve finding the zeros of the order 
parameter (see Section 0) . The density-density correla- 
tor is defined as 



S(r',OHI^M)| 2 hMr + r',i + 0| 2 } 

-(|V;(r,i)| 2 )(hMr + r',t + i')| 2 ) 



(14) 



This correlator and its decay in time is most revealing in 
k space if one is interested in the structure of the system. 
Therefore its Fourier transform is measured, normalized 
by the average density of the order parameter and for 
easier readability divided by its high temperature limit: 



C(M) 



5(k,t) 



|2\2 



lim 



5(k,0) 



(15) 



To compute this quantity the concept of the Fourier 
transform on a plane has to be adapted to the curved 
two dimensional space of the surface of a sphere. The 
Fourier transform is an expansion of the correlation func- 
tion in the complete orthonormal set of functions solving 
the free wave equation, which is in the plane the continu- 
ous set of functions e ik r . The equivalent set of functions 
in a spherical geometry is the discrete set of normalized 
spherical harmonics, Yj m (r). To a value of I corresponds 
k = l/R. Because the system is isotropic, the correlator 
in k space depends only on the magnitude of k, i.e. only 
on / and not on m. So, for better averaging we always 
calculate the correlator for all m and average over the 
different m. 



S(l/R,t) 



21 + 



1 1 

— y 

+i z - 



d 2 rY l m {v)S{v,t). (16) 



Now we substitute for S(l/R, t) from Eq. (|l4|), expand 
tp in LLL eigenfunctions and take the spatial integral in- 
side the thermal average and the summation over the 
lowest Landau levels. The high temperature limit can 
easily be worked out analytically |ll| and the correlation 
function in Eq. ^ can be written in terms of thermal 
averages of the coefficients: 
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FIG. 1. Structure factor for olt = —8, squares represent 
D=0 (pure limit), triangles D=9, and circles D=25. 



C(l/R,t-t') = 



1 



En=0 V n V n) 



4tt(7V ~l)\(N + 1 + 1)1 

(iV!)2(2Z + l) 

/ / min(JV,iV— m) 

E ( E <+m(*K(*)J& 
in— — I \ n — m ax ( , — m ) 

xmn(N,N—m) \ 
n'-max(0,-m) / 



where I™ n is defined in Eq. [l^. Here (...) c sig- 
nifies the connected part of the correlator, i.e. 
(v*(t)v*(t')v k (t) Vl (t')) c = (v*{t) Vl {t'))(v*(t')v k {t)). 
Without disorder the non-connected part of the density- 
density correlator, denoted C nc (k,t), is zero for all ly^O, 
which is equivalent to translational invariance, {\ip(r, t)\ 2 } 
constant. However, in the disordered case there are per- 
manent correlations due to the structure imposed on the 
vortices by the quenched disorder. This is intuitively 
obvious: a permanent local potential will favor specific 
vortex constellations and thus result in permanent cor- 
relations. Therefore (\ip(r,t)\ 2 ) is a function of r and 
C nc (k,t) is not zero. The infinite time limit of the full 
correlator is equal to the non-connected part, for t' — >oo 
(|V(r,t)|^(r + r',i-M')| 2 ) = (|V>(r, t)\ 2 ) (|^(r', t')\ 2 ). 
We measure C(k,t)+C nc (k,t) directly. To retrieve the 
connected part, C nc (k,t) is subtracted. 



A. The structure factor 

The density-density correlator C(k,t) for t' = is the 
structure factor of the system. For the case of no dis- 
order it has been investigated in detail by Dodgson and 
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Moore |{L5| . We give a short summary of their results, 
which we will extend by adding disorder. The structure 
factor shows peaks (which sharpen as the temperature is 
lowered) near the inverse lattice vectors of the triangu- 
lar vortex lattice. A Lorentzian always gives an excellent 
fit for the peak near the first reciprocal lattice vector. 
The associated correlation length £d, (D for density), is 
proportional to the inverse width of this peak at half 
its maximum, (denoted by 5" 1 ), and found to vary as 
C-D oc \a T \l m - 

Fig. |] shows the structure factor at the same effec- 
tive temperature for no, medium and very strong disor- 
der. Disorder flattens the peaks of the structure factor 
and makes it look rather like that of a clean system at a 
higher temperature. Like for the clean case we can fit the 
first peak to a Lorentzian. However, the region around 
the peak where we get a good fit is narrower than in 
the clean case. Due to this limited fitting regime, the 
limited system size and, especially for strong disorder, 
insufficient averaging over disorder, the errors are rather 
large. 



3.5^ 

3 
2.5 

2 
1.5 

1 



2.4 r 

2.2 - 
2 

A, 

1.8 - 

1.6 - 

1.4 - 

-. 1.2 ; 



\\\\ 



•10 



-8 



-7 



0C T 



FIG. 2. The inverse width of the first peak in the struc- 
ture factor depending on —ar with linear fits to the region 
—9 < qt < —6. Filled circles, squares, triangles, empty cir- 
cles and squares correspond to D=0 (pure limit), D=0.25, 1, 
4 and 9 respectively. For D=0 the system size is N=200, for 
D^O, N=72. The inset is reproduced from reference H. It 
shows equivalent results of the parquet graph approximation 
with linear fits for D=0, 0.1, 0.2 and 0.3. 

Fig.^ shows the inverse peak width at half maximum, 
<5 -1 oc £d, from these fits for different degrees of disor- 
der. At high temperatures the presence of disorder does 
not affect the correlation length £d much. With decreas- 



ing temperature the growth of the correlation length is 
reduced. This effect is intuitively obvious because any 
hexagonal order is now opposed not only by thermal 
fluctuations but also by the pinning of vortices to ran- 
domly distributed potential minima. A linear growth of 
(5 _1 with ctT and a decrease of the gradient \d(5~ 1 )/dceT\ 
with increasing disorder has been predicted using a non- 
perturbative, parquet graph approximation to the two 
dimensional LLL system 0. In the inset of Fig. || a plot 
of 5~ 1 (cit) from Ref. jjj is reproduced. We can now com- 
pare d(S~ 1 ) I dax from the parquet graph approximation 
with our simulation results. The absolute values of the 
gradients in the simulation are roughly a factor 2 larger 
than in Ref. 0. However, the relative decrease of the 
gradient due to disorder agree well for Yeo and Moore's 
and our results. If we refer to the linear fits in Fig. || for 
the case of an increase of disorder from D=0 to D=0.5, 
the change in our simulation is only a factor 1.15 larger. 

ln(C) 
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FIG. 3. Typical relaxation behavior for small k (k=0.5, cir- 
cles) and for larger k (k=2.5) without disorder (squares) and 
with disorder (triangles) 



B. Relaxation time scales 

To discuss the relaxation time scales of the density- 
density correlator and their dependence on the wave vec- 
tor, the effective temperature and the strength of disor- 
der, we need to define a relaxation time ro(fc, qt, D). 
Examples of the relaxation behavior of C(k,t) with and 
without disorder are shown in Fig. [|. Without disorder 
ln(C(fc,t)) shows an almost perfectly linear decay with 
time. Only for very small times is the decay faster than 
in a linear fit. In the case of disorder however, the de- 
cay of ln(C(fc, t)) is linear only for later times and has a 
smaller gradient at earlier times. In this case we define 
td from a fit only to the regime of truly linear exponen- 
tial decay. We fit according to C(k,t) cx exp(— i/ro) to 
extract the relaxation time scale tjj. This can be done 



5 



with great accuracy for larger k. However, with and with- 
out disorder, the error in this fit grows large for small k. 
For k < 1.5 we have C(k,t) <C 1. The relaxation times 
are rather small in this regime and C(k 1 t) decays very 
quickly to zero. Very near zero, however, noise dominates 
the data, so that in the worst cases we have to restrict 
our fits to the first two or three points of the curve. 

Having defined td from the exponential decay of the 
density-density correlator, it is an interesting observation 
that the time scales which arise from the exponential de- 
cay of the current correlator (Vxj(r, i)Vxj(r', t 1 )) which 
we measured when determining the transverse conduc- 
tivity ||^] are exactly the same. This suggests these time 
scales are important for all the dynamical properties of 
the system, as is bourne out by our results on the self- 
diffusion coefficient in Section 
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FIG. 4. Logarithmic Relaxation time scales of the decay 
of the density-density correlator in reciprocal space. Circles, 
squares, triangles pointing up and down represent D=0 (pure 
limit) and ar = -12,-10,-8 and -5 respectively. For filled sym- 
bols system size N=200, for empty symbols N=224. Dia- 
monds represent D=4, ar =-8 and N=72. 

We find that, in the vortex liquid regime which we are 
investigating, T£>(fc, ax, D) reflects, like the structure fac- 
tor, the hexagonal lattice structure of the ground state. 
Fig. H shows a logarithmic plot of the k dependence of the 
relaxation times for different temperatures and degrees of 
disorder. Peaks can clearly be seen near the first, second 
and third reciprocal lattice vector of a triangular lattice 
in k space, which lie at k values of 2.694, 4.665 and 5.387 
p3| . Note that for the clean system at ar < —10 the 
second and third peak are clearly resolved, which is not 
the case for the structure factor in the same temperature 
regime. For all measured temperatures, the longest relax- 
ation time occurs at the first peak, i.e. k w G = 2.694. If 



we introduce disorder to the system, the peaks near the 
reciprocal lattice vectors flatten and, like the structure 
factor, the relaxation times look very similar to the ones 
in a clean sample at higher temperature. The error in td 
increases with disorder for all values of k. Like in the case 
of the structure factor, this arises from insufficient aver- 
aging over disorder. The ar dependence of the longest 
time time scales T]j(k « G) can be seen in a logarithmic 
plot in Fig. U for different degrees of disorder. For the 
case without disorder the same quantity is also plotted 
in Fig. ^ over the whole temperature range for which 
measurements were taken. 
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FIG. 5. Logarithm of the largest relaxation time scale ver- 
sus — ar- Filled circles, filled squares, triangles, empty circles, 
and empty squares correspond to D=0, 0.25, 1, 4, and 9 re- 
spectively. For D=0 (pure limit) the system size is N=200, 
for D/0, N=72. 

Without disorder we see an almost perfectly linear ex- 
ponential increase with \olt\-, i- e - t d(G) cx exp(C|aT|)- 
This indicates an activated relaxation process with an 
activation energy Ea/ksT cx ln(r£>) cx \(Xt\- With dis- 
order a general decrease of the relaxation times can be 
observed. The effect of disorder is very small at high tem- 
peratures but grows stronger with decreasing ar- With 
weak disorder the relation \ti(td) cx \olt\ is, at least ap- 
proximately, still valid, and for strong disorder holds also 
within the range of the increasingly large error. Further 
discussion of these results is left to Section VI . 



V. SELF-DIFFUSION 

We have measured the self-diffusion coefficient of the 
vortex positions R. In a liquid the square of the distance 
between initial position of a vortex and its position after 
time t, denoted by S(t), increases linearly in time because 
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of self-diffusion. The self-diffusion coefficient D s is then 
defined from the following relation: 

S{t) = ((R(t) - R(0)) 2 ) = 2D S t. (18) 

In order to compute the self-diffusion coefficient we 
have to monitor the vortex positions on the sphere of 
over not too short a time interval. How the vortex posi- 
tions depend on the LLL coefficients can easily be seen if 
we make a coordinate transformation in Eq. (|7]). We use 
the stereoscopic projection of the sphere into the com- 
plex plane, given by z = x + iy = cos(4>) tan(0/2) + 
isin(0) tan(0/2). Now the order parameter expansion in 
Eq. (Q) can be written as 

O N 

^ )= (^w + i)^ ^" (19) 

The vortices are where ip(z) = 0, i.e. the roots of the Nth 
order polynomial with coefficients v m A m . For a typical 
system size of N=200, the 100th coefficient in these poly- 
nomials is about 60 orders of magnitude larger than the 
first and the last, which makes the numerical root find- 
ing nontrivial. With standard routines to find zeros of 
polynomials failing, we succeeded by searching for zeros 
on a very dense set of trial points using a Laguerre al- 
gorithm (see e.g. Ref. |p6| ). However, we cannot always 
find all vortices on the sphere. Vortices very close to the 
south pole, which is projected to infinity by the stereo- 
scopic projection, can be inaccessible for the numerical 
routine. The number of these inaccessible vortices in- 
creases with system size. With our method, for system 
sizes up to N=250 and for the nearly uniform distribu- 
tion of vortices typically found, the maximum number of 
inaccessible roots is one, and for that one we know that 
it is very near the south pole. 




t/10 3 At 

FIG. 6. <(R(i) - R(0)) 2 ) against t for a T = -7 (stars), 
QfT = —8 (dots), qt = —9 (squares), «t = — 10 (triangles 
pointing up), olt = —11 (triangles pointing down), olt = —12 
(circles) 



Once the roots are found, we have to identify the same 
vortex at different times in order to compute its self- 
diffusion. This was done by simply recording the vortex 
positions in short time intervals At and making a one-to- 
one mapping of the positions at t into the set at t + At 
by pairing the vortices which are the smallest distance 
apart. Only very rarely these mappings were not one-to- 
one, but they could aways be made so by pairing one new 
position with its second or third nearest old position. 



We have recorded the self-diffusion for a system of 200 
vortices and effective temperatures o.t=-7,..,-12 for 13 
(c(t=-7) to 22 (ay=-12) times the density correlation 
relaxation time Td- A plot of S(t) as defined in Eq. ( |l§| ) 
can be seen in Fig. |[ The dependence is linear except for 
very small times t. We have extracted D s from fits to the 
linear regime according to Eq. ([l8]). We are interested 
in how the diffusion time scale tsd oc 1/D s relates to the 
relaxation time scale of the density correlations td . Fig. 

shows both time scales plotted logarithmically against 
oit- The data suggests the same linear exponential de- 
pendence t cx exp(C|aT|) holds for both time scales. If 
we extract the gradient C from the linear fits in Fig. |^ 
for the data for D s and for to for the same system size, 
they agree within 7%. This agreement suggests that self- 
diffusion as well as relaxation processes in the system are 
determined by the same mechanism, with an activation 
energy that grows linearly with — ax- 

ln(l/Db), 
In (x)pv 




FIG. 7. Dominant time scales versus — qt for D=0 (pure 
limit). Empty squares, triangles and crosses: Logarithm 
of the largest relaxation time scale for system sizes of 
N=72 (2 < -a T < 8), 200 (2 < -a T < 12), and 224 
(8 < —otT < 12) respectively. The dotted line is a linear fit for 
N=200. Filled circles: Logarithm of the inverse self-diffusion 
coefficient for N=200 with linear fit (solid line) 







7 



VI. ACTIVATED DYNAMICS 

We now put together our results regarding length and 
time scales in the thin film. Firstly, the length scale of 
translational order £d oc (5 _1 has been measured 15 and 
calculated in the pure limit from the structure factor. 
It increases as £d oc —olt at low temperatures. We found 
that this is at least approximately true in a disordered 
system, a result in good agreement with analytical re- 
sults jjj. Secondly, there is one dominant time scale r 
which determines relaxation and self-diffusion times. It 
can be described as r oc exp(— Colt)- Again this behav- 
ior is not qualitatively changed by disorder. These two 
results together yield 



t oc exp(F£ D /l m ). 



(20) 



Activated dynamics is said to occur when the time 
scale varies exponentially with some power of the cor- 
relation length. This kind of dynamics, which is found to 
be valid for example in spin glasses |3(J is generally de- 
scribed by ln(r) oc ifi . Here L is the linear domain size, 
in our case given by ^d- The barrier-height exponent ip 
describes how the activation energy of the dominant re- 
laxation process, E a , depends on the linear domain size. 
In our system ip=l and E a oc ksT^o- The proportion- 
ality constant F depends on the strength of disorder D. 
In Fig. ||we have plotted ln(r) versus (S^ 1 ) for different 
disorder strengths. The poor quality of our data does 
only allow a quantitative description of F{D). However, 
we can fit the data in the moderate to low temperature 
regime for weak disorder according to Eq. (^) as shown 
in the inset of Fig. |[ A clear decrease of the gradient of 
the linear fits (and hence F) with increasing disorder D 
is visible. 



A. The relaxation mechanism 

Fig. H shows snap shots of the vortex dynamics at 
ax = — 11. The vortex motion is shown over half a self- 
diffusion time tsd(&t)- After this time the average dis- 
placement of a vortex is l m , which is (v / 3/47r) 1 / 2 w 0.37 
times the nearest neighbor distance in the triangular lat- 
tice ground state. The vortex positions are shown for 
24 time steps At between initial and final position. To 
cut out small random moves and thus make the overall 
displacement of a vortex more clearly visible, the posi- 
tions in the picture at time t are an average of the vortex 
positions at times t — At, t and t + At. The different 
grey shades of the vortex positions indicate their coordi- 
nation number. The coordination numbers are not based 
on the time averaged, but on the real positions. The 
nearest-neighbor bonds at initial and final positions are 
also shown. 



8r 
ln(x) 



1.5 2 2.5 , XT 
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FIG. 8. The data from Figs. || and ^ plotted as lnr versus 
<5 _1 . Filled circles, filled squares, triangles, empty circles, and 
empty squares correspond to D=0, 0.25, 1, 4, and 9 respec- 
tively. The inset shows the data for the temperature regime 
—9 < olt < —6 and weak disorder (D=0, 0.25, 1) with linear 
fits. 



To identify pairs of nearest neighbors and calculate 
the coordination numbers we used a "greedy" triangula- 
tion algorithm. It uses the following definition of nearest 
neighbors in non-crystal context: Two vortices are near- 
est neighbors if their connection line does not cross with 
any other belonging to a pair a shorter distance apart. 
To apply this definition to a given configuration, the vor- 
tices are paired in all possible combinations and all pairs 
ordered according to distance. Then it is decided succes- 
sively, starting from the pair with the shortest distance, 
if the vortices of a pair qualify as nearest neighbors, in 
which case they are connected with a bond. A pair is not 
connected with a nearest neighbor bond if its connection 
line crosses over one of the already existing nearest neigh- 
bor bonds. If this is not the case, the vortices in the pair 
are nearest neighbors and a bond is drawn. 

We want to point out the main features about the vor- 
tex dynamics that can be learned from snap shots like the 
one in Fig. |[ At first sight they look quite confusing, 
no easy patterns being distinguishable even after the vor- 
tex paths have been smoothed out. The real relaxation 
mechanisms are not nicely isolated processes in an other- 
wise cryStallirre 'environment. We do for example not see 
"braids" as suggested in reference |l5| , which are formed 
by isolated motion of a ring of vortices. Where motion 
occurs, it does extend over a region of the plane. We also 
see that the vortices in the regions of motion change their 
coordination number very often. Topological defects such 
as dislocations, equivalent to a pair of one fivefold and 
one sevenfold coordinated vortex, occur in large numbers 
and their arrangements may be very transient, especially 
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at higher temperatures. The kind of relaxation process 
we see takes place in a more or less strongly correlated 
liquid and is not well described in terms of isolated topo- 
logical defects in an elastic environment. 



ation times to be larger in larger systems, where these 
defects are less concentrated. Fig. ^ however shows ex- 
cellent agreement of td for systems different in size by 
more than a factor 3. 



0.6 



0.4 



0.2 



-0.2 



-0.4 



-0.6 




-0.6 -0.4 -0.2 0.2 0.4 0. 



FIG. 9. Snap shots of the vortex motion in stereographic 
projection (r=0 is the north pole, r=l is the equator) at 
ar = —11. The nearest-neighbor bonds are in darker grey 
for the starting points and in lighter grey for the end points. 
The vortex positions are shaded in lighter grey if their coordi- 
nation number is <5, in medium grey if it is 6 and in black if 
it is >7. Arrows point out opposite parallel motion of vortex 
chains . 



However, at a second look recurring themes in the vor- 
tex motion become visible. The kind of structure which 
is easiest to spot is a chain-like motion of vortices. The 
process underlying this chain motion can be identified as 
a tilt in the orientation of the nearest-neighbor grid in 
a strongly correlated region of the vortex liquid. In Fig. 
^ such a tilt is easiest to spot as a small clockwise ro- 
tation about the mid point. However, a rotation of the 
underlying grid is not caused by a solid body rotation 
of all the vortices in the region. For a large correlated 
region, a solid body rotation would mean that vortices 
far from the center have to travel a long way, which does 
not happen. Instead each vortex moves in a way to reach 
a position in the new grid which is not too far from its 
initial position. In Fig. [| this effect can be seen in the 
form of an opposite parallel motion of chains of vortices 
indicated with arrows. 

We are confident that we see a genuine property of the 
2D liquid and not a relaxation process favored by the 
fact that on the surface of a sphere there is a permanent 
presence of 12 5-coordinated centers (for detailed discus- 
sion refer to jl^]). If these disclinations were crucial for 
the relaxation mechanism, one should expect the relax- 



B. Activation energy scaling 

We observe that the typical length of moving chains 
increases with decreasing ay. This is easy to understand 
because the range of hexatic order £d oc —olt determines 
the size of a tilt region. This pattern of fewer, large events 
of motion at low temperature and smaller, uncorrelated 
events at higher temperature is confirmed by our self- 
diffusion measurements, where it can be deduced from 
the statistical error in the self-diffusion coefficient; if the 
simulation has been running for the same time in units 
of relaxation times, t/r(o;x)=const., at different ax, the 
vortices in the systems have diffused by the same amount 
but the error is clearly smaller at higher temperatures. 
The range of static correlations in the liquid determines 
the size of a tilt region and therefore the energy barrier 
height of the relaxation mechanism. To explain why the 
activation energy E a depends linearly on £d, we suggest 
the following scaling argument. 

Consider a region of the vortex liquid with linear size 
£d- We want to estimate the energy it costs to tilt the 
core of the region, while the edges relax, as can be seen 
in Fig.^, in order to keep best possible hexatic order with 
the static surroundings. Over a distance £d hexatic cor- 
relations from the middle of the region to the edges are 
only just noticeable. Therefore the core of the region can- 
not rotate freely. The angles which allow some hexatic 
alignment with the surroundings are discrete in steps of 
2tt/N, where N is the number of vortices on the edge 
on the region. The potential energies of the N differ- 
ent orientations with respect to the minimum energy ori- 
entation will have a probability distribution P{E) with 
mean E. A region of size £d is typically in a low energy 
state just locked to a given orientation by its surround- 
ings. The locking energy, which is the lowest potential 
difference E to one of the N other orientations, denoted 
^min, is of order fc^T. By rotating, the system will jump 
from an initial low energy state to higher ones. From 
there, the rotated region and its surroundings relax to 
a new low energy state, which is is the process respon- 
sible for destruction of correlations. The intermediate 
"barrier" state after the rotation is one of the N possible 
orientations, and its energy is a random energy from the 
distribution P(E). We can write the typical activation 
energy as E a = E. From above we know that the typ- 
ical lowest energy in a given sample E m i n rs fc^T. In 
order to find E we can apply the well-known result that 
the minimum value of E which occurs in a large random 
sample of size N is on average E m i n oc E/N . This result 
is shown to be valid for any finite probability distribution 
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P(E), E > with mean E in Ref. §1. We now know 



that E n 



E oc NE„ 



NksT. N is essentially the 



region's perimeter and therefore N oc £d, which yields 
-E^o of C-D the scaling behavior observed in our mea- 
surements. 



VII. CONCLUSIONS 
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We have analyzed the equilibrium dynamics of a thin 
film superconductor in the vortex liquid phase. We have 
shown that there is one main relaxation time scale which 
appears in the decay of density correlations of the order 
parameter as well as in the self-diffusion of the vortices. 
This time scale shows a linear exponential dependence on 
clt at all accessible temperatures. We have also described 
the effects of quenched random disorder on the density 
correlations and their decay. Disorder induces permanent 
correlations and reduces the range of the non-permanent 
correlations and their relaxation times. The density cor- 
relations in a disordered system look very similar to those 
in a system without disorder at a higher temperature. 
This is in good agreement with non-perturbative analyt- 
ical results 0. We find that both the pure and the dis- 
ordered system have activated dynamics. The main time 
scale r scales with the range of translational order £d like 
r oc exp(F £_n/lm), where F depends on the strength of 
disorder. This indicates an activated relaxation process 
with an activation energy of the order E a oc fc^T £i)/Z m . 
We have identified this process as a tilt in the hexagonal 
orientation of the vortex liquid over a region of linear ex- 
tent £d. We have found no divergence of time scales or 
other features attributable to a phase transition at any 
finite temperature. Our work does not confirm the ex- 
istence of a vortex lattice phase at finite temperature, 
nor do we see a vortex glass or Bragg glass phase in the 
disordered system. 

Experimental results on time scales from ac trans- 
port measurement spectra in weakly disordered very films 
should provide information on the relaxation time scales 
found in our numerical work. Such experiments have 
been performed, see for example Refs. |l8) and [|r|. How- 
ever, the samples in thin film experiments are are often 
two dimensional only in the sense that their characteris- 
tic bending length of vortices is larger than the sample 
thickness, but not thin enough to obey 2D LLL scaling. 
Therefore they do not compare directly to our results. 
Experiments similar to the ones from Refs. M] and JTs| ] 
but in thinner films would be directly related to our re- 
sults. Experimental observation of self-diffusion would be 
far more difficult, if not impossible. Decoration experi- 
ments, which use snap shots of the vortex liquid, have 
made detailed studies of static properties possible (see 
e.g. pa] p9|). However, the analysis of the vortex dy- 
namics with this technique remains difficult [E8[ . 
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